MesEpiSubset_2023-11-15

This script is a focused analysis of mesenchymal to epithelial signaling in the BEFM project data, re-run on a deeper standardized sampling (5000 cells per sample), looking specifically at BEF12, BEF14, BEF15, BEFM1 and BEFM2 (per Allie’s request 2023-11-08.)

We will isolate and then look at Mesenchymal Influence, the Epithelial Microenvironment, and Mesencymal-Epithelial Signaling. These three indepednet lenses should allow us to clearly understand what is going on.

This script is a working draft of what might evolve to become a useful standardized framework to look at the full dataset, later. This workflow is experimental for now and is in active revision in collaboration with Allie.

Setup Packages

Here we load our required packages for this analysis.

# Set WD
setwd("/Users/msbr/Library/CloudStorage/GoogleDrive-michasam.raredon@yale.edu/.shortcut-targets-by-id/1VLqBlyzO-Qad5O2kbwXkBRh_1cQho39t/Engineered Lung Paper/Global_Connectomics")

# Require packages
require(Seurat)
require(NICHES)
require(ggplot2)
require(cowplot)
require(dplyr)
require(knitr)

Load NICHES data from Part 1

These data have not been filtered, per Part 2 decisions. This seems to be OK for this project, and generally preferable, because some of our measurements, particularly from BCL5, are low-information and we do not want to lose these data points. So, we load from the last save-point, which contains data from a (maximally) deep sampling of the same number of cells from each sample.

load("~/Library/CloudStorage/GoogleDrive-michasam.raredon@yale.edu/.shortcut-targets-by-id/1VLqBlyzO-Qad5O2kbwXkBRh_1cQho39t/Engineered Lung Paper/Global_Connectomics/global.connectomics.2023-11-11.Robj")

Subset out the just specific samples of interest, per Allie’s request:

Here we limit all of our downstream analysis to information that came exclusively from the following 5 samples of interest, which will allow us to nicely isolate the effect of macrophages on mesenchymal to epithelial signaling.

samples.of.interest <- c('BEF12','BEF14','BEF15','BEFM1','BEFM2')
# we will use this set of samples for downstream calculations in this workflow

Define colors

This is a temporary color palette, and we should consder updating for publication.

load("~/Library/CloudStorage/GoogleDrive-michasam.raredon@yale.edu/.shortcut-targets-by-id/1VLqBlyzO-Qad5O2kbwXkBRh_1cQho39t/Engineered Lung Paper/Allie's Object Filtration/Final Objects/all.color_palettes_2.R")

## For now:
color_pals$type_colors <- c('#A40606','#9CFFFA','#B0DB43','#9C528B','#2F6690',
                                '#946846','#F1C40F','green','#0F0326','#E65F5C','#14591D','#726DA8',
                                'yellow','purple','blue','red','orange','darkgrey','magenta')
names(color_pals$type_colors) <- unique(global.connectomics$CellToCell$alra$SendingType)

Isolate Mesenchymal Influence

This is just the Mesenchymal Cell-To-System data, which incorporate information about both cell-signaling character and cell population presence / absence / representation.

temp <- global.connectomics$CellToSystem$alra # let's just look at the imputed data, for now, for ease
Idents(temp) <- temp$class
table(Idents(temp))

## 
##  Epithelium Endothelium  Mesenchyme      Immune 
##       32888        8164       31702        5071

mes.inf <- subset(temp, idents = 'Mesenchyme')
Idents(mes.inf) <- mes.inf$orig.ident
table(Idents(mes.inf))

## 
##  BEF1  BEF2  BEF3 BEF12 BEF14 BEF15 BEFM1 BEFM2 BEFM4 BEFM5 BEFM6  FB13  FB14 
##  2514   742  1520  4497  4037  4066   986  3250    85    82   143  4878  4902

mes.inf <- subset(mes.inf,idents = c(samples.of.interest))
table(Idents(mes.inf))

## 
## BEF12 BEF14 BEF15 BEFM1 BEFM2 
##  4497  4037  4066   986  3250

Isolate Epithelial Microenvironment

This is just the Epithelial System-To-Cell data, which incorporates information about all signaling landing on Epithelium, regardless of source. It changes, like above, both with changing cellular phenotype and with altered cell population representation.

temp <- global.connectomics$SystemToCell$alra
Idents(temp) <- temp$class
table(Idents(temp))

## 
##  Epithelium Endothelium  Mesenchyme      Immune 
##       32888        8164       31702        5071

epi.mic <- subset(temp, idents = 'Epithelium')
Idents(epi.mic) <- epi.mic$orig.ident
table(Idents(epi.mic))

## 
## BC1P3 BC1P6  BCL5 BCEC2  BEF1  BEF2  BEF3 BEF12 BEF14 BEF15 BEFM1 BEFM2 BEFM4 
##  5000  2887  5000  4967   751   451   632   350   806   781  1502   788   752 
## BEFM5 BEFM6 
##  3875  4346

epi.mic <- subset(epi.mic,idents = c(samples.of.interest))
table(Idents(epi.mic))

## 
## BEF12 BEF14 BEF15 BEFM1 BEFM2 
##   350   806   781  1502   788

Isolate Mesenchymal to Epithelial Signaling

This just looks at the Mesenchymal-To-Epithelium Cell-To-Cell communication. It predominantly captures signal due to altered cellular character. It does not appropriately capture altered cell-signaling cues due to altered cell population representation, but it MAY capture signal present in this cross DUE to the addition of unseen cells outside of this cross (i.e., it could change with/without the presence of macrophages, even though there will be no data FROM macrophages in this object.)

temp <- global.connectomics$CellToCell$alra
Idents(temp) <- temp$class.Joint
table(Idents(temp))

## 
##   Epithelium - Epithelium  Epithelium - Endothelium  Endothelium - Epithelium 
##                     87508                     12249                     12249 
## Endothelium - Endothelium   Epithelium - Mesenchyme  Endothelium - Mesenchyme 
##                     26980                     17089                      8196 
##   Mesenchyme - Epithelium  Mesenchyme - Endothelium   Mesenchyme - Mesenchyme 
##                     17089                      8196                     55356 
##       Epithelium - Immune      Endothelium - Immune       Mesenchyme - Immune 
##                      4343                      3040                      2050 
##       Immune - Epithelium      Immune - Endothelium       Immune - Mesenchyme 
##                      4343                      3040                      2050 
##           Immune - Immune 
##                      9763

mes.epi <- subset(temp,idents = 'Mesenchyme - Epithelium')
Idents(mes.epi) <- mes.epi$orig.ident
table(Idents(mes.epi))

## 
##  BEF1  BEF2  BEF3 BEF12 BEF14 BEF15 BEFM1 BEFM2 BEFM4 BEFM5 BEFM6 
##  2212  1071  1376  1386  2578  2540  2511  2134   360   335   586

mes.epi <- subset(mes.epi,idents = c(samples.of.interest))
table(Idents(mes.epi))

## 
## BEF12 BEF14 BEF15 BEFM1 BEFM2 
##  1386  2578  2540  2511  2134

Ok, now we have isolated our data of interest. Next, let’s visualize each data set to get a first sense of ‘where’ there might be signal in the data that we are interested in exploring more deeply:

Mesenchymal Influence

mes.inf <- ScaleData(mes.inf)
mes.inf <- FindVariableFeatures(mes.inf)
mes.inf <- RunPCA(mes.inf,npcs = 100)
ElbowPlot(mes.inf,ndims = 100)

mes.inf <- RunUMAP(mes.inf,dims = 1:25)

p1 <- DimPlot(mes.inf,group.by = 'Dataset2',cols=color_pals$dataset_colors,raster=F,shuffle = T)
p2 <- DimPlot(mes.inf,group.by = 'orig.ident',cols = color_pals$sample_colors,raster=F,shuffle = T,label=T)
p3 <- DimPlot(mes.inf,group.by = 'class',cols = color_pals$class_colors,raster=F,shuffle = T)
p4 <- DimPlot(mes.inf,group.by = 'SendingType',cols = color_pals$type_colors,raster=F,shuffle = T)

plot_grid(p1,p2,p3,p4)

pdf(file='MesenchymeInfluence_Subset.UMAPS.pdf',width=10,height=8)
print(plot_grid(p1,p2,p3,p4))
dev.off()

## quartz_off_screen 
##                 2

I think this look pretty good. Significant tissue-tissue separation, which may be a consequence of real biology or may be batch, we do not know yet. We could try integrating these later if we would like, not sure if it would help us or obscure interesting patterns. Let’s avoid for the moment, but can revisit this decision.

Note that the above plot shows how fibroblast influence is different between samples, but considers all signals, not just those landing on epithelial cells. A portion of the variance shown here is likely due to changes in signaling landing on epitehlial cells. We do not know yet.

Epithelial Microenvironment

epi.mic <- ScaleData(epi.mic)
epi.mic <- FindVariableFeatures(epi.mic)
epi.mic <- RunPCA(epi.mic,npcs = 100)
ElbowPlot(epi.mic,ndims = 100)

epi.mic <- RunUMAP(epi.mic,dims = 1:25)

p1 <- DimPlot(epi.mic,group.by = 'Dataset2',cols=color_pals$dataset_colors,raster=F,shuffle = T)
p2 <- DimPlot(epi.mic,group.by = 'orig.ident',cols = color_pals$sample_colors,raster=F,shuffle = T,label=T)
p3 <- DimPlot(epi.mic,group.by = 'class',cols = color_pals$class_colors,raster=F,shuffle = T)
p4 <- DimPlot(epi.mic,group.by = 'ReceivingType',cols = color_pals$type_colors,raster=F,shuffle = T)

plot_grid(p1,p2,p3,p4)

pdf(file='EpithelialMicroenviornment_Subset.UMAPS.pdf',width=10,height=8)
print(plot_grid(p1,p2,p3,p4))
dev.off()

## quartz_off_screen 
##                 2

This is actually pretty intresting and I hven’t noticed these patterns before. The tissues group more together here. The similarity between Mono + Tri_E vs. Co + Tri_L + Quad is interesting and I suspect relevant.

BEFM1 is clearly different from the rest, suggesting a relatively unique epithelial microenvironment. This aligns with our biological intuition, given the histology and Allie’s analysis.

Mesenchymal-Epithelial Communication

mes.epi <- ScaleData(mes.epi)
mes.epi <- FindVariableFeatures(mes.epi)
mes.epi <- RunPCA(mes.epi,npcs = 100)
ElbowPlot(mes.epi,ndims = 100)

mes.epi <- RunUMAP(mes.epi,dims = 1:20)

p1 <- DimPlot(mes.epi,group.by = 'Dataset2.Sending',cols=color_pals$dataset_colors,raster=F,shuffle = T)
p2 <- DimPlot(mes.epi,group.by = 'orig.ident',cols = color_pals$sample_colors,raster=F,shuffle = T)
p3 <- DimPlot(mes.epi,group.by = 'class.Sending',cols = color_pals$class_colors,raster=F,shuffle = T)
p4 <- DimPlot(mes.epi,group.by = 'class.Receiving',cols = color_pals$class_colors,raster=F,shuffle = T)
p5 <- DimPlot(mes.epi,group.by = 'SendingType',cols = color_pals$type_colors,raster=F,shuffle = T)
p6 <- DimPlot(mes.epi,group.by = 'ReceivingType',cols = color_pals$type_colors,raster=F,shuffle = T)

plot_grid(p1,p2,p3,p4,p5,p6)

pdf(file='MesenchymalToEpitheliumCommunication_Subset.UMAPS.pdf',width=20,height=10)
print(plot_grid(p1,p2,p3,p4,p5,p6))
dev.off()

## quartz_off_screen 
##                 2

We can also cluster these data, which is rapdily becoming my favorite thing to do:

mes.epi <- FindNeighbors(mes.epi,dims = 1:20)
mes.epi <- FindClusters(mes.epi,res=0.1)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 11149
## Number of edges: 369367
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9313
## Number of communities: 4
## Elapsed time: 1 seconds

p1 <- DimPlot(mes.epi,group.by = 'seurat_clusters')

# Take a quick look at how metadata slots compare
table(mes.epi$seurat_clusters,mes.epi$Dataset2.Sending)

##    
##     Quad_E Tri_L
##   0   1380  5340
##   1   2308    42
##   2    876  1063
##   3     81    59

table(mes.epi$seurat_clusters,mes.epi$orig.ident) # We probably want to take a look at this in a plot

##    
##     BEF12 BEF14 BEF15 BEFM1 BEFM2
##   0  1030  2230  2080    63  1317
##   1     8    13    21  2298    10
##   2   344   299   420   121   755
##   3     4    36    19    29    52

table(mes.epi$SendingType,mes.epi$orig.ident) # Note: No pericyte signaling captured in BEFM1...why? Good thing / Bad thing? Is this insightful/useful information?

##                        
##                         BEF12 BEF14 BEF15 BEFM1 BEFM2
##   Aberrant_Fibroblast     350   806   620   165   788
##   Alveolar_Fibroblast     350   806   781  1502   788
##   Cycling_Mes             350   651   730   769   214
##   Pericytes               116    91    50     0   176
##   Remodeling_Fibroblast   220   224   359    75   168

# This measures how different archetypes are distributed over each tissue sample
to.plot <- data.frame(table(mes.epi$seurat_clusters,mes.epi$orig.ident))
p2 <- ggplot(to.plot,
       aes(x=Var2,y=Freq,fill=Var1,group=Var1))+
  #geom_bar(stat='identity',position = 'dodge')+
  theme_classic()+
  geom_density(stat = 'identity',alpha=0.5)+
  ggtitle('NICHES Archetype Distribution Over Sample')

# This measures how different sending cell populations are distributed over each tissue sample
p3 <- DimPlot(mes.epi,group.by = 'SendingType',cols = color_pals$type_colors,raster=F,shuffle = T)
to.plot <- data.frame(table(mes.epi$SendingType,mes.epi$orig.ident))
# Order the factors for plotting
temp <- to.plot %>% group_by(Var1) %>% summarise(Frequency = sum(Freq))
temp <- arrange(temp,desc(Frequency))
to.plot$Var1 <- factor(to.plot$Var1,levels = temp$Var1)
p4 <- ggplot(to.plot,
       aes(x=Var2,y=Freq,fill=Var1,group=Var1))+
  theme_classic()+
  geom_density(stat = 'identity',alpha=0.5,)+
  scale_fill_manual(values=color_pals$type_colors)+
  ggtitle('SendingType Distribution Over Sample')

# This measures how different receiving cell populations are distributed over each tissue sample
p5 <- DimPlot(mes.epi,group.by = 'ReceivingType',cols = color_pals$type_colors,raster=F,shuffle = T)
to.plot <- data.frame(table(mes.epi$ReceivingType,mes.epi$orig.ident))
# Order the factors for plotting
temp <- to.plot %>% group_by(Var1) %>% summarise(Frequency = sum(Freq))
temp <- arrange(temp,desc(Frequency))
to.plot$Var1 <- factor(to.plot$Var1,levels = temp$Var1)
p6 <- ggplot(to.plot,
       aes(x=Var2,y=Freq,fill=Var1,group=Var1))+
  theme_classic()+
  geom_density(stat = 'identity',alpha=0.5,)+
  scale_fill_manual(values=color_pals$type_colors)+
  ggtitle('ReceivingType Distribution Over Sample')

plot_grid(p1,p3,p5,p2,p4,p6)

pdf(file='MesenchymalToEpitheliumCommunication_Subset_Clustered_Distributions.pdf',width=20,height=10)
plot_grid(p1,p3,p5,p2,p4,p6)
dev.off()

## quartz_off_screen 
##                 2

I think this is a cool new plotting technique. What do you think?

This is very cool! Looks like there is a clear effect on Mesenchymal-Epithelial communication due to the addition of Macrophages. Let’s dig into this specific signal a bit by running a statistical test comparing these two categories.

Couple things here. 1. There is a clear association between alveolar fibroblasts with Quad culture, that is already known I think from your analysis but it is cool to see here. 2. Seems like aberrant fibroblasts go down with the addition of macrophages, which again is known from your analysis but i also think is cool. 3. There is a small but possibly real increase in BEFM1 of signaling landing on BASCs. What does this mean / how do we interpret? 4. The BEFM1 archetype seems to be deeply correlated with Alveolar_Fibroblast development. 5. I think that signaling archetype 0 is associated with the Aberrant_Fibroblast phenotype. Let’s check:

# This measures how different sending cell populations are distributed over each signaling archetype
to.plot <- data.frame(table(mes.epi$SendingType,mes.epi$seurat_clusters))
# Order the factors for plotting
temp <- to.plot %>% group_by(Var1) %>% summarise(Frequency = sum(Freq))
temp <- arrange(temp,desc(Frequency))
to.plot$Var1 <- factor(to.plot$Var1,levels = temp$Var1)

plot <- ggplot(to.plot,
       aes(x=Var2,y=Freq,fill=Var1,group=Var1))+
  theme_classic()+
  geom_density(stat = 'identity',alpha=0.5,)+
  scale_fill_manual(values=color_pals$type_colors)+
  ggtitle('SendingType Distribution Over Archetype')+
  xlab('Signaling Archetype')
plot_grid(p4,plot,nrow=1)

Huh. Not what I was expecting. I think instead, what we are seeing is that Archetype 1 preferentially does not include Aberrant_Fibroblasts (Aberrant_Fibroblasts do not generally communicate with epithelium using the mechanisms present in signaling Archetype 1.) We already know that Archetype 1 is associated with BEFM1. Very cool!!

There are a couple of ways that we can parse out specific signals that we might find to be important biologically. We can leverage any meta-data slot(s) that we want to perform these calculations We could directly compare Tri Culture vs. Quad Culture within this data set, as follows:

Compare MesEpi with and without macrophages:

First, we want a new meta data category that annotates Tri culture vs. Quad culture. we do this as follows:

mes.epi$Dataset4 <- NA
Idents(mes.epi) <- mes.epi$Dataset2.Sending 
mes.epi <- RenameIdents(mes.epi,
                        'Tri_E'='Tri',
                        'Tri_L'='Tri',
                        'Quad_E'='Quad',
                        'Quad_L'='Quad')
mes.epi$Dataset4 <- Idents(mes.epi)
table(mes.epi$Dataset4)

## 
##  Tri Quad 
## 6504 4645

Now, we can run a marker test comparing these two categories specifically:

mark.mes.epi <- FindAllMarkers(mes.epi,only.pos = T,min.pct = 0.1,logfc.threshold = 0.01)
mark.mes.epi$ratio <- mark.mes.epi$pct.1/mark.mes.epi$pct.2
mark.mes.epi$power <- mark.mes.epi$ratio*mark.mes.epi$avg_log2FC
kable(mark.mes.epi %>% group_by(cluster) %>% top_n(50,ratio) %>% arrange(desc(ratio)))

p_val	avg_log2FC	pct.1	pct.2	p_val_adj	cluster	gene	ratio	power
0	0.0406835	0.255	0.000	0	Quad	Wnt7a—Lrp6	Inf	Inf
0	0.0359223	0.117	0.001	0	Quad	Sct—Sctr	117.000000	4.2029065
0	0.3779125	0.558	0.014	0	Quad	Adipoq—Adipor1	39.857143	15.0625137
0	0.0848352	0.209	0.006	0	Quad	Lamb3—Itga2	34.833333	2.9550926
0	0.2268704	0.326	0.012	0	Quad	Adipoq—Adipor2	27.166667	6.1633118
0	0.0766549	0.225	0.010	0	Quad	Lamc2—Itga2	22.500000	1.7247360
0	0.0433286	0.194	0.010	0	Quad	Lama1—Itga2	19.400000	0.8405755
0	0.0257216	0.167	0.009	0	Quad	Tnfsf12—Tnfrsf25	18.555556	0.4772778
0	0.0279528	0.230	0.013	0	Quad	Csf2—Il3ra	17.692308	0.4945493
0	0.0598572	0.192	0.012	0	Quad	Cdh1—Egfr	16.000000	0.9577156
0	0.2540416	0.374	0.024	0	Quad	Lamb3—Itga3	15.583333	3.9588143
0	0.0269311	0.139	0.009	0	Quad	Csf2—Csf2ra	15.444444	0.4159356
0	0.0452318	0.126	0.009	0	Tri	Rtn4—Tnfrsf19	14.000000	0.6332454
0	0.0724622	0.206	0.015	0	Quad	S100a9—Tlr4	13.733333	0.9951475
0	0.0425042	0.143	0.012	0	Quad	Cdh1—Igf1r	11.916667	0.5065085
0	0.1358085	0.357	0.030	0	Quad	Lama1—Itga3	11.900000	1.6161206
0	0.2117413	0.373	0.033	0	Tri	Gdf10—Acvr1b	11.303030	2.3933186
0	0.2256925	0.396	0.036	0	Quad	Lamc2—Itga3	11.000000	2.4826174
0	0.0253900	0.198	0.018	0	Quad	Fgf1—Egfr	11.000000	0.2792896
0	0.0202926	0.110	0.010	0	Quad	Tcn2—Lrp2	11.000000	0.2232181
0	0.1194167	0.247	0.023	0	Quad	Lamb3—Itga6	10.739130	1.2824314
0	0.1925922	0.349	0.033	0	Quad	Lamb3—Itgb4	10.575758	2.0368086
0	0.0657472	0.116	0.011	0	Quad	Col5a3—Sdc3	10.545454	0.6933345
0	0.0538167	0.175	0.017	0	Quad	Wnt7b—Lrp5	10.294118	0.5539958
0	0.3899451	0.457	0.045	0	Quad	Lamb3—Itgb1	10.155556	3.9601089
0	0.0262222	0.101	0.010	0	Quad	Mdk—Lrp2	10.100000	0.2648437
0	0.1282287	0.371	0.038	0	Quad	Lama1—Itgb4	9.763158	1.2519167
0	0.0466438	0.107	0.011	0	Quad	Apoe—Lrp2	9.727273	0.4537166
0	0.1443302	0.435	0.046	0	Quad	Lama1—Sdc4	9.456522	1.3648620
0	0.0668396	0.217	0.023	0	Quad	Cdh1—Lrp5	9.434783	0.6306173
0	0.0565876	0.113	0.012	0	Quad	Mdk—Sdc3	9.416667	0.5328662
0	0.0128408	0.129	0.014	0	Quad	Fgf1—Fgfr3	9.214286	0.1183191
0	0.0857056	0.284	0.031	0	Quad	Angptl4—Tie1	9.161290	0.7851740
0	0.0901725	0.279	0.031	0	Quad	Lama1—Itga6	9.000000	0.8115523
0	0.1540670	0.301	0.034	0	Quad	Csf2—Itgb1	8.852941	1.3639458
0	0.2538694	0.494	0.058	0	Quad	Lama1—Itgb1	8.517241	2.1622666
0	0.0529281	0.110	0.013	0	Quad	Wnt5a—Mcam	8.461538	0.4478533
0	0.1192592	0.253	0.030	0	Quad	Cdh1—Ptprf	8.433333	1.0057523
0	0.0723921	0.124	0.015	0	Quad	Col5a1—Sdc3	8.266667	0.5984410
0	0.1745612	0.179	0.022	0	Tri	Slit1—Sdc1	8.136364	1.4202936
0	0.0319312	0.160	0.020	0	Quad	Psen1—Notch4	8.000000	0.2554493
0	0.1649228	0.367	0.046	0	Quad	Lamc2—Itgb4	7.978261	1.3157968
0	0.1016286	0.261	0.033	0	Quad	Lamc2—Itga6	7.909091	0.8037901
0	0.3426242	0.479	0.061	0	Quad	Lamc2—Itgb1	7.852459	2.6904429
0	0.0823475	0.101	0.013	0	Quad	Anxa1—Dysf	7.769231	0.6397771
0	0.0524192	0.246	0.032	0	Quad	Adam15—Itga5	7.687500	0.4029727
0	0.0837317	0.160	0.021	0	Quad	Tgfb1—Cxcr4	7.619048	0.6379560
0	0.0227634	0.107	0.015	0	Quad	Adm—Calcrl	7.133333	0.1623789
0	0.0461006	0.260	0.037	0	Quad	Fgf1—Cd44	7.027027	0.3239505
0	0.2001768	0.219	0.032	0	Quad	Mmp9—Ephb2	6.843750	1.3699600
0	0.0230529	0.116	0.017	0	Quad	Sema3f—Nrp1	6.823529	0.1573019
0	0.0513693	0.199	0.030	0	Quad	Cdh1—Erbb3	6.633333	0.3407498
0	0.0170333	0.131	0.021	0	Quad	Lama1—Gpc1	6.238095	0.1062552
0	0.0322468	0.110	0.026	0	Tri	Thbs1—Itga2b	4.230769	0.1364290
0	0.1465410	0.139	0.033	0	Tri	Fn1—Itga2b	4.212121	0.6172484
0	0.1324566	0.139	0.033	0	Tri	Calr—Itga2b	4.212121	0.5579234
0	0.0750097	0.168	0.043	0	Tri	Hdc—Hrh4	3.906977	0.2930610
0	0.0461880	0.206	0.054	0	Tri	Slit3—Robo2	3.814815	0.1761988
0	0.0516654	0.247	0.066	0	Tri	Slit2—Robo2	3.742424	0.1933539
0	0.3357056	0.228	0.067	0	Tri	Fn1—Nt5e	3.402985	1.1424013
0	0.2020212	0.345	0.104	0	Tri	Col14a1—Cd44	3.317308	0.6701666
0	0.0358461	0.111	0.034	0	Tri	Pigf—Flt1	3.264706	0.1170270
0	0.0320742	0.127	0.039	0	Tri	Ngf—Ngfr	3.256410	0.1044467
0	0.1303773	0.189	0.059	0	Tri	B2m—Cd247	3.203390	0.4176493
0	0.0578606	0.139	0.044	0	Tri	Gnai2—Ptpru	3.159091	0.1827870
0	0.0516398	0.147	0.051	0	Tri	Tnc—Nt5e	2.882353	0.1488441
0	0.1464161	0.335	0.127	0	Tri	Bdnf—Ngfr	2.637795	0.3862157
0	0.2108339	0.129	0.050	0	Tri	Fn1—Itga8	2.580000	0.5439515
0	0.1786970	0.361	0.140	0	Tri	App—Ngfr	2.578571	0.4607830
0	0.1963006	0.360	0.140	0	Tri	Rtn4—Ngfr	2.571429	0.5047729
0	0.6967932	0.781	0.318	0	Tri	Thbs2—Itgb1	2.455975	1.7113067
0	0.5210727	0.732	0.301	0	Tri	Thbs2—Cd47	2.431894	1.2671935
0	0.0356187	0.245	0.104	0	Tri	Mst1—Mst1r	2.355769	0.0839094
0	0.0710921	0.254	0.110	0	Tri	Efemp1—Egfr	2.309091	0.1641582
0	0.0902954	0.129	0.056	0	Tri	A2m—Lrp1	2.303571	0.2080018
0	0.4032791	0.414	0.183	0	Tri	Thbs2—Itga6	2.262295	0.9123364
0	0.1230577	0.166	0.079	0	Tri	Dlk1—Notch1	2.101266	0.2585769
0	0.0882803	0.182	0.088	0	Tri	Dlk1—Notch2	2.068182	0.1825798
0	0.0172538	0.153	0.075	0	Tri	Bmp8b—Bmpr2	2.040000	0.0351978
0	0.1340253	0.370	0.182	0	Tri	Icam4—Itgb1	2.032967	0.2724691
0	0.2346354	0.370	0.189	0	Tri	Fn1—Itgb7	1.957672	0.4593391
0	0.0166589	0.174	0.096	0	Tri	Bmp8b—Acvr2a	1.812500	0.0301943
0	0.0221861	0.175	0.099	0	Tri	Bmp8b—Bmpr1a	1.767677	0.0392179
0	0.2945381	0.512	0.290	0	Tri	Sfrp1—Fzd6	1.765517	0.5200120
0	0.0117863	0.152	0.087	0	Tri	Bmp8b—Acvr1	1.747126	0.0205921
0	0.1396410	0.255	0.147	0	Tri	Efemp2—Lingo1	1.734694	0.2422344
0	0.0260303	0.109	0.063	0	Tri	Apoe—Lrp8	1.730159	0.0450365
0	0.7234671	0.818	0.473	0	Tri	Ptn—Sdc1	1.729387	1.2511545
0	0.0966098	0.254	0.147	0	Tri	Rtn4—Lingo1	1.727891	0.1669312
0	0.6257268	0.762	0.442	0	Tri	Ptn—Plxnb2	1.723982	1.0787417
0	0.0619612	0.517	0.302	0	Tri	Rgma—Bmpr2	1.711920	0.1060726
0	0.1020057	0.241	0.141	0	Tri	Vcam1—Itgb1	1.709220	0.1743501
0	0.1064315	0.545	0.327	0	Tri	Bmp5—Bmpr2	1.666667	0.1773859
0	0.0132838	0.125	0.075	0	Tri	Dkk2—Kremen2	1.666667	0.0221396
0	0.0640236	0.124	0.078	0	Tri	Fn1—Itgb8	1.589744	0.1017811
0	0.0234753	0.124	0.078	0	Tri	Col4a1—Itgb8	1.589744	0.0373197
0	0.1158296	0.335	0.213	0	Tri	Gnai2—S1pr5	1.572770	0.1821733
0	0.0591788	0.260	0.166	0	Tri	Gnai2—Adcy9	1.566265	0.0926896
0	0.2718472	0.655	0.424	0	Tri	Fbln1—Itgb1	1.544811	0.4199527
0	0.1321240	0.598	0.389	0	Tri	Bmp5—Bmpr1a	1.537275	0.2031109

Not bad! Some very interesting signals in there. Let’s remove signaling mediated by integrins, which are less interesting for this project.

MOI <- mark.mes.epi[-grep('Itg',mark.mes.epi$gene),]
MOI <- MOI %>% group_by(cluster) %>% top_n(50,ratio) %>% arrange(desc(ratio))
kable(MOI)

p_val	avg_log2FC	pct.1	pct.2	p_val_adj	cluster	gene	ratio	power
0	0.0406835	0.255	0.000	0	Quad	Wnt7a—Lrp6	Inf	Inf
0	0.0359223	0.117	0.001	0	Quad	Sct—Sctr	117.000000	4.2029065
0	0.3779125	0.558	0.014	0	Quad	Adipoq—Adipor1	39.857143	15.0625137
0	0.2268704	0.326	0.012	0	Quad	Adipoq—Adipor2	27.166667	6.1633118
0	0.0257216	0.167	0.009	0	Quad	Tnfsf12—Tnfrsf25	18.555556	0.4772778
0	0.0279528	0.230	0.013	0	Quad	Csf2—Il3ra	17.692308	0.4945493
0	0.0598572	0.192	0.012	0	Quad	Cdh1—Egfr	16.000000	0.9577156
0	0.0269311	0.139	0.009	0	Quad	Csf2—Csf2ra	15.444444	0.4159356
0	0.0452318	0.126	0.009	0	Tri	Rtn4—Tnfrsf19	14.000000	0.6332454
0	0.0724622	0.206	0.015	0	Quad	S100a9—Tlr4	13.733333	0.9951475
0	0.0425042	0.143	0.012	0	Quad	Cdh1—Igf1r	11.916667	0.5065085
0	0.2117413	0.373	0.033	0	Tri	Gdf10—Acvr1b	11.303030	2.3933186
0	0.0253900	0.198	0.018	0	Quad	Fgf1—Egfr	11.000000	0.2792896
0	0.0202926	0.110	0.010	0	Quad	Tcn2—Lrp2	11.000000	0.2232181
0	0.0657472	0.116	0.011	0	Quad	Col5a3—Sdc3	10.545454	0.6933345
0	0.0538167	0.175	0.017	0	Quad	Wnt7b—Lrp5	10.294118	0.5539958
0	0.0262222	0.101	0.010	0	Quad	Mdk—Lrp2	10.100000	0.2648437
0	0.0466438	0.107	0.011	0	Quad	Apoe—Lrp2	9.727273	0.4537166
0	0.1443302	0.435	0.046	0	Quad	Lama1—Sdc4	9.456522	1.3648620
0	0.0668396	0.217	0.023	0	Quad	Cdh1—Lrp5	9.434783	0.6306173
0	0.0565876	0.113	0.012	0	Quad	Mdk—Sdc3	9.416667	0.5328662
0	0.0128408	0.129	0.014	0	Quad	Fgf1—Fgfr3	9.214286	0.1183191
0	0.0857056	0.284	0.031	0	Quad	Angptl4—Tie1	9.161290	0.7851740
0	0.0529281	0.110	0.013	0	Quad	Wnt5a—Mcam	8.461538	0.4478533
0	0.1192592	0.253	0.030	0	Quad	Cdh1—Ptprf	8.433333	1.0057523
0	0.0723921	0.124	0.015	0	Quad	Col5a1—Sdc3	8.266667	0.5984410
0	0.1745612	0.179	0.022	0	Tri	Slit1—Sdc1	8.136364	1.4202936
0	0.0319312	0.160	0.020	0	Quad	Psen1—Notch4	8.000000	0.2554493
0	0.0823475	0.101	0.013	0	Quad	Anxa1—Dysf	7.769231	0.6397771
0	0.0837317	0.160	0.021	0	Quad	Tgfb1—Cxcr4	7.619048	0.6379560
0	0.0227634	0.107	0.015	0	Quad	Adm—Calcrl	7.133333	0.1623789
0	0.0461006	0.260	0.037	0	Quad	Fgf1—Cd44	7.027027	0.3239505
0	0.2001768	0.219	0.032	0	Quad	Mmp9—Ephb2	6.843750	1.3699600
0	0.0230529	0.116	0.017	0	Quad	Sema3f—Nrp1	6.823529	0.1573019
0	0.0513693	0.199	0.030	0	Quad	Cdh1—Erbb3	6.633333	0.3407498
0	0.0170333	0.131	0.021	0	Quad	Lama1—Gpc1	6.238095	0.1062552
0	0.0349009	0.111	0.018	0	Quad	Serpine1—Lrp2	6.166667	0.2152221
0	0.0141602	0.117	0.019	0	Quad	Bmp3—Acvr1	6.157895	0.0871973
0	0.0295976	0.134	0.022	0	Quad	Bmp3—Bmpr1a	6.090909	0.1802765
0	0.0846650	0.224	0.038	0	Quad	S100a8—Tlr4	5.894737	0.4990782
0	0.0242458	0.127	0.022	0	Quad	Bmp3—Acvr2a	5.772727	0.1399646
0	0.0221524	0.114	0.020	0	Quad	Lrpap1—Lrp2	5.700000	0.1262684
0	0.0537176	0.296	0.053	0	Quad	Angpt1—Tie1	5.584906	0.3000077
0	0.0880642	0.111	0.020	0	Quad	Ptn—Sdc3	5.550000	0.4887564
0	0.0157701	0.100	0.019	0	Quad	Bmp3—Bmpr2	5.263158	0.0830005
0	0.0202024	0.146	0.028	0	Quad	Angptl2—Tie1	5.214286	0.1053409
0	0.0980915	0.191	0.040	0	Quad	Mmp9—Lrp1	4.775000	0.4683869
0	0.1495824	0.343	0.072	0	Quad	Tgfb1—Eng	4.763889	0.7125938
0	0.0145269	0.128	0.027	0	Quad	Fgf1—Fgfr2	4.740741	0.0688684
0	0.2366882	0.315	0.067	0	Quad	Mmp9—Cd44	4.701492	1.1127878
0	0.0628092	0.233	0.053	0	Quad	Lin7c—Abca1	4.396226	0.2761235
0	0.2203756	0.366	0.085	0	Quad	Hbegf—Egfr	4.305882	0.9489113
0	0.0714967	0.262	0.061	0	Quad	Ntn1—Unc5b	4.295082	0.3070844
0	0.0750097	0.168	0.043	0	Tri	Hdc—Hrh4	3.906977	0.2930610
0	0.0461880	0.206	0.054	0	Tri	Slit3—Robo2	3.814815	0.1761988
0	0.0516654	0.247	0.066	0	Tri	Slit2—Robo2	3.742424	0.1933539
0	0.3357056	0.228	0.067	0	Tri	Fn1—Nt5e	3.402985	1.1424013
0	0.2020212	0.345	0.104	0	Tri	Col14a1—Cd44	3.317308	0.6701666
0	0.0358461	0.111	0.034	0	Tri	Pigf—Flt1	3.264706	0.1170270
0	0.0320742	0.127	0.039	0	Tri	Ngf—Ngfr	3.256410	0.1044467
0	0.1303773	0.189	0.059	0	Tri	B2m—Cd247	3.203390	0.4176493
0	0.0578606	0.139	0.044	0	Tri	Gnai2—Ptpru	3.159091	0.1827870
0	0.0516398	0.147	0.051	0	Tri	Tnc—Nt5e	2.882353	0.1488441
0	0.1464161	0.335	0.127	0	Tri	Bdnf—Ngfr	2.637795	0.3862157
0	0.1786970	0.361	0.140	0	Tri	App—Ngfr	2.578571	0.4607830
0	0.1963006	0.360	0.140	0	Tri	Rtn4—Ngfr	2.571429	0.5047729
0	0.5210727	0.732	0.301	0	Tri	Thbs2—Cd47	2.431894	1.2671935
0	0.0356187	0.245	0.104	0	Tri	Mst1—Mst1r	2.355769	0.0839094
0	0.0710921	0.254	0.110	0	Tri	Efemp1—Egfr	2.309091	0.1641582
0	0.0902954	0.129	0.056	0	Tri	A2m—Lrp1	2.303571	0.2080018
0	0.1230577	0.166	0.079	0	Tri	Dlk1—Notch1	2.101266	0.2585769
0	0.0882803	0.182	0.088	0	Tri	Dlk1—Notch2	2.068182	0.1825798
0	0.0172538	0.153	0.075	0	Tri	Bmp8b—Bmpr2	2.040000	0.0351978
0	0.0166589	0.174	0.096	0	Tri	Bmp8b—Acvr2a	1.812500	0.0301943
0	0.0221861	0.175	0.099	0	Tri	Bmp8b—Bmpr1a	1.767677	0.0392179
0	0.2945381	0.512	0.290	0	Tri	Sfrp1—Fzd6	1.765517	0.5200120
0	0.0117863	0.152	0.087	0	Tri	Bmp8b—Acvr1	1.747126	0.0205921
0	0.1396410	0.255	0.147	0	Tri	Efemp2—Lingo1	1.734694	0.2422344
0	0.0260303	0.109	0.063	0	Tri	Apoe—Lrp8	1.730159	0.0450365
0	0.7234671	0.818	0.473	0	Tri	Ptn—Sdc1	1.729387	1.2511545
0	0.0966098	0.254	0.147	0	Tri	Rtn4—Lingo1	1.727891	0.1669312
0	0.6257268	0.762	0.442	0	Tri	Ptn—Plxnb2	1.723982	1.0787417
0	0.0619612	0.517	0.302	0	Tri	Rgma—Bmpr2	1.711920	0.1060726
0	0.1064315	0.545	0.327	0	Tri	Bmp5—Bmpr2	1.666667	0.1773859
0	0.0132838	0.125	0.075	0	Tri	Dkk2—Kremen2	1.666667	0.0221396
0	0.1158296	0.335	0.213	0	Tri	Gnai2—S1pr5	1.572770	0.1821733
0	0.0591788	0.260	0.166	0	Tri	Gnai2—Adcy9	1.566265	0.0926896
0	0.1321240	0.598	0.389	0	Tri	Bmp5—Bmpr1a	1.537275	0.2031109
0	0.3591713	0.585	0.384	0	Tri	Gpc3—Cd81	1.523438	0.5471751
0	0.1389964	0.364	0.241	0	Tri	Fgf10—Fgfr2	1.510373	0.2099365
0	0.0275825	0.221	0.147	0	Tri	Gnai2—Cxcr2	1.503401	0.0414676
0	0.0779625	0.199	0.134	0	Tri	Dcn—Met	1.485075	0.1157801
0	0.0802552	0.601	0.405	0	Tri	Bmp5—Acvr2a	1.483951	0.1190947
0	0.0402334	0.293	0.200	0	Tri	Bmp5—Acvr2b	1.465000	0.0589419
0	0.1508341	0.361	0.248	0	Tri	Apoe—Ldlr	1.455645	0.2195609
0	0.0754248	0.523	0.360	0	Tri	Bmp5—Acvr1	1.452778	0.1095755
0	0.0171711	0.193	0.134	0	Tri	Cntf—Il6r	1.440298	0.0247316
0	0.0694165	0.218	0.153	0	Tri	B2m—Cd3g	1.424837	0.0989072
0	0.0560992	0.190	0.136	0	Tri	Cxcl12—Ackr3	1.397059	0.0783739
0	0.1617107	0.433	0.314	0	Tri	Cthrc1—Fzd6	1.378981	0.2229959

So there is definitely a huge activation of WNT family signaling here. Wnt7a, Wnt7b, via Lrp5/6. Rspo1-Lgr4 (I know the ligand well, but not the receptor!) Also mild upregulation of Fgf1/2/7 signaling via Fgfr3. Also a big upregulation of Hbegf-Egfr. The tri-culture elevation of Dlk1-Notch1 is interesting, not sure what it means. Nrg1-Erbb3 also upregulated in Quad, and Erbb3 is an powerful epithelial stem-cell / growth activator. The Csf2-Csf2ra upregulated in Quad might be relevant to macrophage biology. Efna4—Epha7 also upregulated, it is a spatial guidance molecule implicated in tissue morphogenesis, as is Sema4b—Dcbld2, which is acting through a receptor I have never seen before. Il11—Il6st also upregulated. Immune recruitment / activation? Not sure what this mechanism does, would need to read.

Let’s make some feature plots of Mes-Epi signals that are upregulated in Quad culture as compared to Tri, like this:

quick.function <- function(mech){
p1 <- FeaturePlot(mes.epi,mech)
p2 <- VlnPlot(mes.epi,mech)
p3 <- VlnPlot(mes.epi,mech,group.by = 'orig.ident')
plot_grid(p1,p2,p3,nrow = 1)
}
mech <- 'Fgf7—Fgfr3'
quick.function(mech)

mech <- 'Wnt7a—Lrp6'
quick.function(mech)

mech <- 'Wnt7b—Lrp5'
quick.function(mech)

mech <- 'Csf2—Csf2ra'
quick.function(mech)

mech <- 'Fgf2—Fgfr3'
quick.function(mech)

mech <- 'Fgf2—Fgfr2'
quick.function(mech)

mech <- 'Hbegf—Egfr'
quick.function(mech)

mech <- 'Il11—Il6st'
quick.function(mech)

mech <- 'Rspo1—Lgr4'
quick.function(mech)

EDITS 11/9/2023, post-subset: I have added a third plot here that looks at the data split by Tissue, so that we can see if any of the effects are sample-specific. So, while these mechanisms are definitely still interesting, I am seeing the same thing you have mentioned, which is that BEFM1 is like “much better” than other quad cultures. Weird. But I guess also makes sense? We already kind of knew this. Not a huge issue, the tissue looks the best too…

Next, let’s look at some Tri culture Mes-Epi markers:

Not sure what Dlk1 is doing or if it is relevant, it is a Notch ligand and it is definitely lowered in BEFM1:

mech <- 'Dlk1—Notch1'
quick.function(mech)

Thbs2 is Themis’ pro-fibrotic, anti-angiogenic matrix molecule that he knocks out, wild that it is here so clearly and is so reduced in BEFM1 and the quad cultures more generally

mech <- 'Thbs2—Cd47'
quick.function(mech)

Sfrp1—Fzd6 is an extracellular binder of WNT ligands that sequesters them outside of cells and rpevents them from binding, it is net anti-WNT. It is DOWNREGULATED in quad culture, suggesting that macrophage addition will activate this pathway and promote epithelial WNT signal transduction

mech <-'Sfrp1—Fzd6'
quick.function(mech)

Cthrc1—Fzd6 is a major matrix molecule implicated in fibrosis (ligand): https://www.nature.com/articles/s41467-020-15647-5 and the receptor is a critical cue for alveolar regeneration by epithelium: https://www.cell.com/cell/pdf/S0092-8674(23)00541-X.pdf

mech <-'Cthrc1—Fzd6'
quick.function(mech)

But I prefer to look at how the signaling archetypes difer from one another. It is more generalizable, and also it leverages the actual data structure to perform the entire analysis, which means that the signals that come out as output mechanisms are much clearer observed in the data itself. It is easier to interpret this way. So, here is a parallel workflow which uses the archetype information as the main data handle:

Differential testing by Archetype

Idents(mes.epi) <- mes.epi$seurat_clusters # This makes the archetype information the main handle
mark.mes.epi <- FindAllMarkers(mes.epi,only.pos = T,min.pct = 0.1,logfc.threshold = 0.01)
mark.mes.epi$ratio <- mark.mes.epi$pct.1/mark.mes.epi$pct.2
mark.mes.epi$power <- mark.mes.epi$ratio*mark.mes.epi$avg_log2FC
kable(mark.mes.epi %>% group_by(cluster) %>% top_n(10,ratio))

p_val	avg_log2FC	pct.1	pct.2	p_val_adj	cluster	gene	ratio	power
0	0.3664661	0.391	0.084	0	0	App—Ngfr	4.654762	1.7058126
0	0.3373300	0.390	0.084	0	0	Rtn4—Ngfr	4.642857	1.5661749
0	0.2120839	0.361	0.078	0	0	Bdnf—Ngfr	4.628205	0.9815679
0	0.0599403	0.278	0.048	0	0	Mst1—Mst1r	5.791667	0.3471544
0	0.1153768	0.278	0.057	0	0	Gnai2—Cxcr2	4.877193	0.5627149
0	0.0899865	0.221	0.047	0	0	Cxcl1—Cxcr2	4.702128	0.4231281
0	0.0399730	0.215	0.049	0	0	Vasp—Cxcr2	4.387755	0.1753916
0	0.0743570	0.150	0.034	0	0	Dkk2—Kremen2	4.411765	0.3280457
0	0.0379441	0.129	0.032	0	0	Ngf—Ngfr	4.031250	0.1529620
0	0.0262909	0.110	0.025	0	0	Vcam1—Itgb7	4.400000	0.1156799
0	0.6116117	0.941	0.054	0	1	Adipoq—Adipor1	17.425926	10.6579010
0	0.4329424	0.651	0.041	0	1	Lamb3—Itga3	15.878049	6.8742809
0	0.3752132	0.541	0.036	0	1	Adipoq—Adipor2	15.027778	5.6386208
0	0.1516485	0.373	0.015	0	1	Lamb3—Itga2	24.866667	3.7709915
0	0.1404703	0.409	0.017	0	1	Lamc2—Itga2	24.058823	3.3795490
0	0.0736566	0.475	0.008	0	1	Wnt7a—Lrp6	59.375000	4.3733602
0	0.0681361	0.220	0.003	0	1	Sct—Sctr	73.333333	4.9966447
0	0.0470222	0.241	0.016	0	1	Csf2—Csf2ra	15.062500	0.7082716
0	0.0469633	0.356	0.023	0	1	Fgf1—Egfr	15.478261	0.7269106
0	0.0543063	0.120	0.007	0	1	Calcb—Ramp1	17.142857	0.9309652
0	0.1040387	0.294	0.024	0	2	Wnt5a—Ror1	12.250000	1.2744745
0	0.0834049	0.272	0.013	0	2	Cntf—Lifr	20.923077	1.7450864
0	0.0631832	0.238	0.022	0	2	Pon2—Htr2a	10.818182	0.6835278
0	0.2816088	0.157	0.003	0	2	Fn1—Col13a1	52.333333	14.7375259
0	0.1553227	0.205	0.017	0	2	Gnai2—Adra2b	12.058824	1.8730091
0	0.1008910	0.140	0.003	0	2	Nid1—Col13a1	46.666667	4.7082451
0	0.0531654	0.137	0.002	0	2	Nid2—Col13a1	68.500000	3.6418282
0	0.0658717	0.141	0.004	0	2	Icam1—Il2ra	35.250000	2.3219768
0	0.0269528	0.135	0.009	0	2	Ctf1—Lifr	15.000000	0.4042924
0	0.1260287	0.110	0.010	0	2	Col1a1—Itga11	11.000000	1.3863162
0	0.6905106	0.907	0.007	0	3	Vegfa—Itga9	129.571429	89.4704476
0	0.6463646	0.921	0.007	0	3	Spp1—Itga9	131.571429	85.0431074
0	0.3968352	0.693	0.005	0	3	Tnc—Itga9	138.600000	55.0013534
0	0.3476353	0.929	0.007	0	3	Vegfc—Itga9	132.714286	46.1361757
0	0.2938716	0.664	0.004	0	3	Adam12—Itga9	166.000000	48.7826932
0	0.0902893	0.343	0.002	0	3	Adam15—Itga9	171.500000	15.4846157
0	0.0677749	0.186	0.001	0	3	Vcam1—Itga9	186.000000	12.6061392
0	0.0207233	0.136	0.001	0	3	Csf2—Itga9	136.000000	2.8183676
0	0.0125385	0.150	0.001	0	3	Wnt7a—Fzd9	150.000000	1.8807717
0	0.0308627	0.129	0.001	0	3	Lamb3—Cd151	129.000000	3.9812839

When we then look at feature plots, things look very nice. We can view trends in the data that correlate with sample, but we didn’t use those labels to perform the analysis. This keeps the analysis indepedent and non-biased.

FeaturePlot(mes.epi,'Wnt7a—Lrp6',max.cutoff = 0.2,split.by = 'orig.ident')

FeaturePlot(mes.epi,'Csf2—Csf2ra',max.cutoff = 0.5,split.by = 'orig.ident',order = T)

FeaturePlot(mes.epi,'Sct—Sctr',max.cutoff = 0.5,split.by = 'orig.ident',order = T) # What is this molecular interaction? It's wild.

FeaturePlot(mes.epi,'Cxcl1—Cxcr2',max.cutoff = 0.9,split.by = 'orig.ident',order=T) # BEFM1 has reduced inflammation ????

FeaturePlot(mes.epi,'Wnt5a—Ror1',max.cutoff = 0.9,split.by = 'orig.ident',order=T)

Reduced Wnt5a in BEFM ???? I think of Wnt5a as a “pro-differentation” cue…could the BEFM1 condition be promoting stemness in the epithelial cells? Is this possible / supported by your analysis?

Downregulated Wnt5a and upregulated Wnt7a … connected to the increased BASC representation in this archetype?

If we want to, we can create a complex heatmap which shows a great deal of information at once, something like as follows:

Let’s make a heatmap that displays all key findings for Mes-Epi Signaling in this dataset

source('/Users/msbr/Library/CloudStorage/GoogleDrive-michasam.raredon@yale.edu/My Drive/Tuft_Sox9_Pneumonectomy_Project/Homeostatic_Single_Cell/CustomHeatmap.R')

mechs.to.annotate <- mark.mes.epi %>% group_by(cluster) %>% top_n(10,ratio)
CustomHeatmap(object = mes.epi,
              data.type = 'CellToCell',
              primary = 'seurat_clusters' ,
              secondary = 'SendingType' ,
              tertiary = 'ReceivingType' ,
              quarternary = 'orig.ident' ,
              primary.cols = NULL,
              secondary.cols = color_pals$type_colors, # Need to be a named list of colors in the right order
              tertiary.cols = color_pals$type_colors,
              quarternary.cols = color_pals$sample_colors,
              features = unique(mark.mes.epi$gene),
              labels = c('Signaling Archetype','Sending Cell Type','Receiving Cell Type','Tissue'),
              selected.row.anotations = mechs.to.annotate$gene,
              selected.label.size = 10)

As a weird sort of cross-project scientific experiment, what happens if I ask to annotate code-inheritied markers of interest from the pneumonectomy project??

mechs.to.annotate <- c('Il17b—Il17rb','Il25—Il17rb',
                                          'Ereg—Erbb3','Areg—Erbb3',
                                          'Sema4d—Plxnb2','Sema6a—Plxna2',#'Sema6a—Plxna4',
                                          'Wnt3a—Fzd6','Wnt7b—Lrp5','Wnt5a—Fzd5','Wnt5a—Fzd6','Wnt3a—Ryk',
                                          'Rspo2—Lgr5','Rspo3—Sdc4','Rspo1—Lrp6','Rspo1—Lrp6','Rspo1—Lgr5',
                                          'Dll1—Notch2','Dll4—Notch2','Dll3—Notch2',
                                          'Sct—Sctr','Dhh—Ptch1',
                                          'Tgfb1—Tgfbr1','Tgfb1—Tgfbr2','Ccl2—Ackr4','Ccl5—Sdc1',#'Ccl5—Sdc4',
                                          'Mdk—Sdc4','Ccl11—Ackr4','Fbln1—Itgb1',
                                          'Il18—Cd48','Tnf—Ltbr','Igf1—Insr','Tnfsf13—Tnfrsf13b')
CustomHeatmap(object = mes.epi,
              data.type = 'CellToCell',
              primary = 'seurat_clusters' ,
              secondary = 'SendingType' ,
              tertiary = 'ReceivingType' ,
              quarternary = 'orig.ident' ,
              primary.cols = NULL,
              secondary.cols = color_pals$type_colors, # Need to be a named list of colors in the right order
              tertiary.cols = color_pals$type_colors,
              quarternary.cols = color_pals$sample_colors,
              features = unique(mark.mes.epi$gene),
              labels = c('Signaling Archetype','Sending Cell Type','Receiving Cell Type','Tissue'),
              selected.row.anotations = mechs.to.annotate,
              selected.label.size = 10)

Fascinating. Seems like there are some markers that overlap and might be biologcally relevant in both cases. Could we learn from this? Could we adopt a similar apporach to compare engineered biology to native? I need to think on this more. Seems powerful to me though.

My overarching summary, in a nutshell:

Tri-culture has major upregualtion of fibrotic pathways, genes, and matrix molecules. Adding macrophages basically removes this signal.

Adding macrophages in Quad-culture promotes significant upregulation of WNT-signaling and FGF-family signaling. We also see the expression of prominent sptial-guidance molecules crucial to local tissue organziation that we know from normal alveolar homeostasis.

I’m beginning to consider if the BEFM1 condition / Quad-Culture condition / addition of macrophages might promote stemness in these cultures. Would this make sense with the rest of our data?